Chapter 8 - Exponential smoothing

R. Hyndman/R. J. Serrano

7/30/2022

Exponential smoothing

Historical perspective

Big idea: control the rate of change

\(\alpha\) controls the flexibility of the level

\(\beta\) controls the flexibility of the trend

\(\gamma\) controls the flexibility of the seasonality

ETS models

Additive ("A") or multiplicative ("M")

None ("N"), additive ("A"), multiplicative ("M"), or damped ("Ad" or "Md").

None ("N"), additive ("A") or multiplicative ("M")

Simple exponential smoothing

Simple methods

Time series \(y_1,y_2,\dots,y_T\).

Simple Exponential Smoothing

Simple Exponential Smoothing

Algerian exports ETS alpha animation

Simple Exponential Smoothing

Iterate to get exponentially weighted moving average form.

Optimising smoothing parameters

Simple Exponential Smoothing

Models and methods

Methods

Models

ETS(A,N,N): SES with additive errors

Forecast error: \(e_t = y_t - \pred{y}{t}{t-1} = y_t - \ell_{t-1}\).

Specify probability distribution for \(e_t\), we assume \(e_t = \varepsilon_t\sim\text{NID}(0,\sigma^2)\).

ETS(A,N,N): SES with additive errors

where \(\varepsilon_t\sim\text{NID}(0,\sigma^2)\).

ETS(M,N,N): SES with multiplicative errors.

ETS(A,N,N): Specifying the model

ETS(y ~ error("A") + trend("N") + season("N"))

By default, an optimal value for \(\alpha\) and \(\ell_0\) is used.

\(\alpha\) can be chosen manually in trend().

trend("N", alpha = 0.5)
trend("N", alpha_range = c(0.2, 0.8))

Example: Algerian Exports

algeria_economy <- global_economy %>%
  filter(Country == "Algeria")

fit <- algeria_economy %>%
  model(ANN = ETS(Exports ~ error("A") + trend("N") + season("N")))
report(fit)
## Series: Exports 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.84 
## 
##   Initial states:
##  l[0]
##  39.5
## 
##   sigma^2:  35.6
## 
##  AIC AICc  BIC 
##  447  447  453

Example: Algerian Exports

components(fit) %>% autoplot()

Example: Algerian Exports

components(fit) %>%
  left_join(fitted(fit), by = c("Country", ".model", "Year"))

Example: Algerian Exports

fit %>%
  forecast(h = 5) %>%
  autoplot(algeria_economy) +
  labs(y = "% of GDP", title = "Exports: Algeria")

Models with trend

Holt’s linear trend

ETS(A,A,N)

Holt’s linear method with additive errors.

Exponential smoothing: trend/slope

Australian population trend beta animation

ETS(M,A,N)

Holt’s linear method with multiplicative errors.

ETS(A,A,N): Specifying the model

ETS(y ~ error("A") + trend("A") + season("N"))

By default, optimal values for \(\beta\) and \(b_0\) are used.

\(\beta\) can be chosen manually in trend().

trend("A", beta = 0.004)
trend("A", beta_range = c(0, 0.1))

Example: Australian population

aus_economy <- global_economy %>% filter(Code == "AUS") %>%
  mutate(Pop = Population / 1e6)
fit <- aus_economy %>%
  model(AAN = ETS(Pop ~ error("A") + trend("A") + season("N")))
report(fit)
## Series: Pop 
## Model: ETS(A,A,N) 
##   Smoothing parameters:
##     alpha = 1 
##     beta  = 0.327 
## 
##   Initial states:
##  l[0]  b[0]
##  10.1 0.222
## 
##   sigma^2:  0.0041
## 
##   AIC  AICc   BIC 
## -77.0 -75.8 -66.7

Example: Australian population

components(fit) %>% autoplot()

Example: Australian population

components(fit) %>%
  left_join(fitted(fit), by = c("Country", ".model", "Year"))

Example: Australian population

fit %>%
  forecast(h = 10) %>%
  autoplot(aus_economy) +
  labs(y = "Millions", title = "Population: Australia")

Damped trend method

Example: Australian population

aus_economy %>%
  model(holt = ETS(Pop ~ error("A") + trend("Ad") + season("N"))) %>%
  forecast(h = 20) %>%
  autoplot(aus_economy)

Example: Australian population

fit <- aus_economy %>%
  filter(Year <= 2010) %>%
  model(
    ses = ETS(Pop ~ error("A") + trend("N") + season("N")),
    holt = ETS(Pop ~ error("A") + trend("A") + season("N")),
    damped = ETS(Pop ~ error("A") + trend("Ad") + season("N"))
  )
tidy(fit)
accuracy(fit)

Example: Australian population

term SES Linear trend Damped trend
\(\alpha\) 1.00 1.00 1.00
\(\beta^*\) 0.30 0.40
\(\phi\) 0.98
NA 0.22 0.25
NA 10.28 10.05 10.04
Training RMSE 0.24 0.06 0.07
Test RMSE 1.63 0.15 0.21
Test MASE 6.18 0.55 0.75
Test MAPE 6.09 0.55 0.74
Test MAE 1.45 0.13 0.18

Models with seasonality

Holt-Winters additive method

Holt and Winters extended Holt’s method to capture seasonality.

Holt-Winters additive method

Exponential smoothing: seasonality

Medicare Australia cost of vaccine scripts seasonality gamma animation

ETS(A,A,A)

Holt-Winters additive method with additive errors.

Holt-Winters multiplicative method

Seasonal variations change in proportion to the level of the series.

ETS(M,A,M)

Holt-Winters multiplicative method with multiplicative errors.

Example: Australian holiday tourism

aus_holidays <- tourism %>%
  filter(Purpose == "Holiday") %>%
  summarise(Trips = sum(Trips))
fit <- aus_holidays %>%
  model(
    additive = ETS(Trips ~ error("A") + trend("A") + season("A")),
    multiplicative = ETS(Trips ~ error("M") + trend("A") + season("M"))
  )
fc <- fit %>% forecast()

Example: Australian holiday tourism

fc %>%
  autoplot(aus_holidays, level = NULL) +
  labs(y = "Thousands", title = "Overnight trips")

Estimated components

components(fit)

Estimated components

Holt-Winters damped method

Often the single most accurate forecasting method for seasonal data:

Holt-Winters with daily data

sth_cross_ped <- pedestrian %>%
  filter(
    Date >= "2016-07-01",
    Sensor == "Southern Cross Station"
  ) %>%
  index_by(Date) %>%
  summarise(Count = sum(Count) / 1000)
sth_cross_ped %>%
  filter(Date <= "2016-07-31") %>%
  model(
    hw = ETS(Count ~ error("M") + trend("Ad") + season("M"))
  ) %>%
  forecast(h = "2 weeks") %>%
  autoplot(sth_cross_ped %>% filter(Date <= "2016-08-14")) +
  labs(
    title = "Daily traffic: Southern Cross",
    y = "Pedestrians ('000)"
  )

Holt-Winters with daily data

Innovations state space models

Exponential smoothing methods

ETS models

Additive error models

Multiplicative error models

Estimating ETS models

Innovations state space models

Let \(\bm{x}_t = (\ell_t, b_t, s_t, s_{t-1}, \dots, s_{t-m+1})\) and \(\varepsilon_t\stackrel{\mbox{\scriptsize iid}}{\sim} \mbox{N}(0,\sigma^2)\).
Additive errors
\(k(x)=1\).\(y_t = \mu_{t} + \varepsilon_t\).
Multiplicative errors
\(k(\bm{x}_{t-1}) = \mu_{t}\).\(y_t = \mu_{t}(1 + \varepsilon_t)\). \(\varepsilon_t = (y_t - \mu_t)/\mu_t\) is relative error.

Innovations state space models

Parameter restrictions

Usual region

Admissible region

Model selection

where \(L\) is the likelihood and \(k\) is the number of parameters initial states estimated in the model.

which is the AIC corrected (for small sample bias).

AIC and cross-validation

Automatic forecasting

From Hyndman et al. (IJF, 2002):

Method performed very well in M3 competition.

Example: National populations

fit <- global_economy %>%
  mutate(Pop = Population / 1e6) %>%
  model(ets = ETS(Pop))
fit

Example: National populations

fit %>%
  forecast(h = 5)

Example: Australian holiday tourism

holidays <- tourism %>%
  filter(Purpose == "Holiday")
fit <- holidays %>% model(ets = ETS(Trips))
fit

Example: Australian holiday tourism

fit %>%
  filter(Region == "Snowy Mountains") %>%
  report()
## Series: Trips 
## Model: ETS(M,N,A) 
##   Smoothing parameters:
##     alpha = 0.157 
##     gamma = 1e-04 
## 
##   Initial states:
##  l[0] s[0] s[-1] s[-2] s[-3]
##   142  -61   131 -42.2 -27.7
## 
##   sigma^2:  0.0388
## 
##  AIC AICc  BIC 
##  852  854  869

Example: Australian holiday tourism

fit %>%
  filter(Region == "Snowy Mountains") %>%
  components(fit)

Example: Australian holiday tourism

fit %>%
  filter(Region == "Snowy Mountains") %>%
  components(fit) %>%
  autoplot()

Example: Australian holiday tourism

fit %>% forecast()

Example: Australian holiday tourism

fit %>% forecast() %>%
  filter(Region == "Snowy Mountains") %>%
  autoplot(holidays) +
  labs(y = "Thousands", title = "Overnight trips")

Residuals

Response residuals

\[\hat{e}_t = y_t - \hat{y}_{t|t-1}\]

Innovation residuals

Additive error model: \[\hat\varepsilon_t = y_t - \hat{y}_{t|t-1}\]

Multiplicative error model: \[\hat\varepsilon_t = \frac{y_t - \hat{y}_{t|t-1}}{\hat{y}_{t|t-1}}\]

Example: Australian holiday tourism

aus_holidays <- tourism %>%
  filter(Purpose == "Holiday") %>%
  summarise(Trips = sum(Trips))
fit <- aus_holidays %>%
  model(ets = ETS(Trips)) %>%
  report()
## Series: Trips 
## Model: ETS(M,N,M) 
##   Smoothing parameters:
##     alpha = 0.358 
##     gamma = 0.000969 
## 
##   Initial states:
##  l[0]  s[0] s[-1] s[-2] s[-3]
##  9667 0.943 0.927 0.968  1.16
## 
##   sigma^2:  0.0022
## 
##  AIC AICc  BIC 
## 1331 1333 1348

Example: Australian holiday tourism

residuals(fit)
residuals(fit, type = "response")

Example: Australian holiday tourism

fit %>%
  augment()

Some unstable models

Exponential smoothing models

Forecasting with exponential smoothing

Forecasting with ETS models

iterate the equations for \(t=T+1,T+2,\dots,T+h\) and set all \(\varepsilon_t=0\) for \(t>T\).

Example: ETS(A,A,N)

\[\begin{align*} y_{T+1} &= \ell_T + b_T + \varepsilon_{T+1}\\ \hat{y}_{T+1|T} & = \ell_{T}+b_{T}\\ y_{T+2} & = \ell_{T+1} + b_{T+1} + \varepsilon_{T+2}\\ & = (\ell_T + b_T + \alpha\varepsilon_{T+1}) + (b_T + \beta \varepsilon_{T+1}) + \varepsilon_{T+2} \\ \hat{y}_{T+2|T} &= \ell_{T}+2b_{T} \end{align*}\] etc.

Example: ETS(M,A,N)

\[\begin{align*} y_{T+1} &= (\ell_T + b_T )(1+ \varepsilon_{T+1})\\ \hat{y}_{T+1|T} & = \ell_{T}+b_{T}.\\ y_{T+2} & = (\ell_{T+1} + b_{T+1})(1 + \varepsilon_{T+2})\\ & = \left\{ (\ell_T + b_T) (1+ \alpha\varepsilon_{T+1}) + \left[b_T + \beta (\ell_T + b_T)\varepsilon_{T+1}\right] \right\} (1 + \varepsilon_{T+2}) \\ \hat{y}_{T+2|T} &= \ell_{T}+2b_{T} \end{align*}\] etc.

Forecasting with ETS models

can only be generated using the models.

Prediction intervals

Example: Corticosteroid drug sales

h02 <- PBS %>%
  filter(ATC2 == "H02") %>%
  summarise(Cost = sum(Cost))
h02 %>% autoplot(Cost)

Example: Corticosteroid drug sales

h02 %>%
  model(ETS(Cost)) %>%
  report()
## Series: Cost 
## Model: ETS(M,Ad,M) 
##   Smoothing parameters:
##     alpha = 0.307 
##     beta  = 0.000101 
##     gamma = 0.000101 
##     phi   = 0.978 
## 
##   Initial states:
##    l[0] b[0]  s[0] s[-1] s[-2] s[-3] s[-4] s[-5] s[-6] s[-7] s[-8] s[-9]
##  417269 8206 0.872 0.826 0.756 0.773 0.687  1.28  1.32  1.18  1.16   1.1
##  s[-10] s[-11]
##    1.05  0.981
## 
##   sigma^2:  0.0046
## 
##  AIC AICc  BIC 
## 5515 5519 5575

Example: Corticosteroid drug sales

h02 %>%
  model(ETS(Cost ~ error("A") + trend("A") + season("A"))) %>%
  report()
## Series: Cost 
## Model: ETS(A,A,A) 
##   Smoothing parameters:
##     alpha = 0.17 
##     beta  = 0.00631 
##     gamma = 0.455 
## 
##   Initial states:
##    l[0] b[0]   s[0]   s[-1]   s[-2]   s[-3]   s[-4]  s[-5]  s[-6]  s[-7]
##  409706 9097 -99075 -136602 -191496 -174531 -241437 210644 244644 145368
##   s[-8] s[-9] s[-10] s[-11]
##  130570 84458  39132 -11674
## 
##   sigma^2:  3.5e+09
## 
##  AIC AICc  BIC 
## 5585 5589 5642

Example: Corticosteroid drug sales

h02 %>%
  model(ETS(Cost)) %>%
  forecast() %>%
  autoplot(h02)

Example: Corticosteroid drug sales

h02 %>%
  model(
    auto = ETS(Cost),
    AAA = ETS(Cost ~ error("A") + trend("A") + season("A"))
  ) %>%
  accuracy()
Model MAE RMSE MAPE MASE RMSSE
auto 38649 51102 4.99 0.638 0.689
AAA 43378 56784 6.05 0.716 0.766